Local density of states and Priedel oscillations in graphene 
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We investigate the local density of states and Friedel oscillation in graphene around a well localized 
impurity in Born approximation. In our analytical calculations Green's function technique has been 
used taking into account both the localized atomic wavefunctions in a tight-binding scheme and the 
corresponding symmetries of the lattice. As a result we obtained long wavelength oscillations in 
the density of electrons with long range behavior proportional to the inverse square of the distance 
from the impurity. These leading oscillations are out of phase on nearby lattice sites (in fact for 
an extended defect they cancel each other within one unit cell), therefore a probe with resolution 
worse than a few unit cells will experience only the next to leading inverse cube decay of density 
oscillations even for a short range scatterer. 
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Ever since the first production of atomically thin car- 
bon filmsi, graphene continues to fascinate physicists of 
almost all walks of life. This simple two dimensional 
system of carbon atoms arranged in a honeycomb lat- 
tice provides us with a number of interesting properties 
due mainly to the massless Dirac nature of the disper- 
sion relation of its electrons^. Undoped graphene has 
its Fermi energy at the tip of the Dirac cones and be- 
haves as a zero gap semiconductor. Applying appropriate 
gate voltage leads to electron or hole pockets and metal- 
lic behavior with Fermi wavenumber kp typically much 
smaller than the size of the Brillouin zone. Exciting po- 
tential applications of graphene include for example car- 
bon based planar electronic circuitry with the possibility 
of electrically reconfigurable wiring 3 -, exploiting the guid- 
ing effect of graphene p-n junctions with negative refrac- 
tive index^. On the theoretical side perhaps the simplest 
quantity to be considered is the change in the local den- 
sity of states (LDOS) and the Friedel oscillation (FO) in 
the excess charge density due to a well localized impu- 
rity. Results of these considerations have implications for 
the LDOS in disordered graphene^, for the interaction 
between adatoms in graphene 6 , or in case of magnetic 
adatomsl for the corresponding RKKY interaction 8 . 

Early theoretical work on the LDOS and the resulting 
FO around an impurity in graphene^^ predicted long 
wavelength (2k f) oscillations in the charge density, but 
with envelope decaying like r -3 at distance r from the 
impurity. This is in contrast to the r~ 2 decay in a de- 
generate nonrelativistic two dimensional Fermi gas, and 
suppressed backscattering of chiral graphene electrons 
residing around the Fermi circle of the Dirac cone was 
offered as an explanation. However, graphene has two 
inequivalent Dirac cones (valleys) in the Brillouin zone, 
and intervalley scattering by the impurity may lead to 
short wavelength oscillations on the order of a few lat- 
tice constants as well. Indeed, a scanning tunneling mi- 
croscopy (STM) studyii of epitaxial graphene revealed 
two different length scales around defects. Subsequently, 
intervalley (or internodal) scattering has been built in 



the theoryi 2 -, and has been shown to be responsible for 
r _1 decay in LDOS and r~ 2 decay in FO. The different 
power laws for intranodal and internodal scattering have 
been observed by Fourier Transform STM (FT-STM) 
experiment^ 3 -, but no detailed investigation of the short 
wavelength oscillations has been performed either exper- 
imentally or theoretically. The first step towards this 
direction has been made by incorporating atomic wave- 
functions instead of just plane waves into the theory^, 
and it has become clear that although the LDOS falls 
off with r , for intranodal scattering only it has oppo- 
site sign on the two sublattices. Therefore upon coarse- 
grainig, for example due to experimental resolution, the 
leading order decay cancels within one unit cell, and one 
is left with the next to leading r~ 2 envelope, and conse- 
quently an r -3 decay of FO. 

The aim of the present paper is to give a detailed anal- 
ysis of the short wavelength oscillations due to internodal 
scattering within the framework of a simple and transpar- 
ent, atomic resolution theory based on the tight-binding 
wavefunctions of electrons on the honeycomb lattice. Our 
main results are the short range spatial pattern of LDOS 
and FO on both sublattices with the corresponding sym- 
metries (see Fig. [1]), and the observation that includ- 
ing internodal scattering still lead to cancellation of the 
leading power law decay of oscillations, but not within 
one unit cell, but within three neighboring ones. This 
means that if for example the STM tip can not resolve 
the six atoms on a hexagon of the honeycomb lattice, 
we encounter again just the next to leading power law 
envelopes for both LDOS and FO. 

We begin with the tight-binding Hamiltonian of carbon 
atoms forming a honeycomb lattice with nearest neighbor 
hopping t (assumed to be real) between the p z orbitals 
tp(r). In momentum space this leads to a 2x2 Hamilto- 
nian matrix spanned by Bloch waves constructed from 
the atomic orbitals on the two sublattices A and B as 
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where /(k) = 1 + e* kai +e lka2 with primitive lattice vec- 
tors ai : 2 = a(±l/2, \/3/2) of length a. Diagonalization 
yields the eigenvalues e±(k) = ±|i||/(k)| defining the 
two bands of pure graphene, and the eigenvectors in the 
tight-binding form 
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where the summation runs over the R lattice vectors and 
6 (It) is the complex phase of /(k) determining the mixing 
of A and B sublattice states, and c = (ax + a 2 )/3 is the 
vector pointing from the A to the B site within the unit 
cell. As is well known, /(k) vanishes at the corners of the 
hexagonal Brillouin zone, leading to the massless Dirac 
spectrum e±(k) ~ ihvpSk , with vp — \fZa\t\j2K the 
Fermi velocity, and 5k the wavenumber measured from 
the corners. Since there are two inequivalent corners, for 
example K = (bi — b 2 )/3 and K' = -K as expressed 
by the primitive reciprocal lattice vectors, there are two 
Dirac cones or valleys in the Brillouin zone. Electron 
spin does not play any role in the forthcoming discussion, 
therefore we suppress spin indices, and all our results will 
refer to one spin orientation. 

The real space representation of the Green's func- 
tion of complex energy variable G(z, r, r') is of central 
importance for our subject, since the LDOS is given 
by its analytic continuation just above the real axis as 
p(s,r) = — 7r _1 ImG(e + id, r, r). For pure graphene, not 
perturbed by impurities, the free Green's function is eval- 
uated using Eq.@ as 
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where f(r) = [y(r), <p(r — c)] is a row vector formed by 
the orbitals of the A and B sites in the unit cell, and the 
Green's matrix is given by 



the proximity of r and r' to a lattice point given by ei- 
ther ta(= R) on the A sublattice, or by r B (= R + c) on 
the B sublattice. On the other hand, the (identical) di- 
agonal elements of the Green's matrix in Eq.Q are even 
functions of R (eg. G AA (z, -R) = G° AA (z,R)), while 
the off diagonal elements transform into each other upon 
reflection (G AB (z,— R) = G BA (z,H)). It is instruc- 
tive to consider the LDOS of the pure system given by 
p (e,r) = Po (£)ER[Mr-R)| 2 + |^(r-R-c)| 2 ]/2, where 
Po(e) is the DOS, which can be approximated around the 
Dirac point by A c \e\/ir(hvF) 2 , A c = \/3a 2 /2 being the 
area of the unit cell. 

Let us consider a substitutional impurity characterized 
by a short range scattering potential energy U(r) = uS(r) 
located at the origin, which is a site on the A sublattice. 
The correction to the LDOS in Born approximation is de- 
termined by AG(z,r, r) = uG°(z, r, 0)G°(z, 0, r), which 
can be expressed due to the localized atomic orbitals as 



u \ip(r-r A )\ 2 G° AA (z,r A )G° AA (z,-r A ), or (5) 
u \ip(r - r B )\ 2 G% A (z,T B - c)G° AB (z, -r B + c) , (6) 

depending on whether r is in the vicinity of an A or a B 
site respectively, and uo = u\ip(0)\ 2 . Clearly the spatial 
pattern is dominated by the density profile of the atomic 
orbital centered on the given lattice site, while the multi- 
plicative factor on that site depends on the corresponding 
Green's matrix elements. It can be proven that the di- 
agonal elements of the Green's matrix like G aa (z,Ta) 
have sixfold rotational symmetry in r A , while the off di- 
agonal elements like G BA (z,r B — c) have only threefold 
rotational symmetry in r B . 

In order to evaluate the relevant Green's matrix ele- 
ments from Eq.Q we observe that the most important 
contributions come from the nodal points of the spec- 
trum, i.e. from around the points K and K' = — K 
in the Brillouin zone. The matrix elements therefore 
consist of two terms each, led by fast oscillations of 
the type exp(iKR) and exp(iK'R), modulated by func- 
tions of slow spatial variation. These latter functions 
can be evaluated by using the linearized /(±K + <5k) = 
— \/3a(±5k x +i5ky) /2 expression around the nodal points 
up to a cutoff k c . This will be a good approximation 
for these slowly varying factors of the Green's matrix 
elements for distances from the impurity much larger, 
and for characteristic spatial variations much longer than 
l/k c . This procedure leads to the following result for the 
diagonal element in Eq.©: 
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The free Green's function in Eq.([3]) is clearly translation 
invariant by a lattice vector only, but due to the well 
localized nature of the atomic wavefunctions it typically 
consists of only one non-negligible term determined by 
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where Hq 1 \z) is the Hankel function, and this functional 
form is valid for Imz > only (which is enough for the 
prescribed analytic continuation) and for \z\ <C hvpk c . 
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Similarly, for the off diagonal element in Eq. (j6j) we obtain 
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where ijj^(z) = ~dH^\z) / dz 1 {} is the angle 
makes with the a; axis (which is parallel to K). For 
G^ B (z, — vb + c) we have the same formula as in Eq.©, 
except that the prefactor with the exponentials is re- 
placed by (_ e -i(Kr B +.9) + e -i(K'r B -#)), We re mind the 

reader that these exponential prefactors in the Green's 
matrix elements describing short wavelength oscillations 
are exact, as opposed to the spatial dependence described 
by the Hankel functions. 

The change of the LDOS in graphene due to the im- 
purity Ap(e,r) = — Tr~ 1 ImAG(e + iS, r,r) is now easily 
calculated using Eqs.([5][5]). Let us note first, that the 
short wavelength spatial pattern in eg. Eq.© is given 
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However, if the impurity is unable to produce intervalley 
scattering (eg. because it is extended and has no large 
wavenumber Fourier components), then plane waves be- 
longing to different valleys will not contribute to the 
above product, and we will only have 1 + 1 = 2 as a 
result. Thus, without intervalley scattering, the factor 
cos 2 (Kr J 4) describing the short wavelength spatial pat- 
tern should be replaced by its average i.e. 1/2. The same 
is true for the factor sin 2 (Kr^ + $) appearing in Eq. © . 

Returning to the change in the LDOS, after analytic 
continuation we obtain 
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if r is in the vicinity of an A site, and 
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if r is in the vicinity of a B site. Here Jo, J\ and Yo, 
are Bessel functions of the first and the second kind. The 
factors describing the short wavelength spatial behavior 
are given by 

c A {v) = \<p(T - r A )\ 2 cos 2 (Kr A ) , and (11) 
c B (r) = ^(r-r B )| 2 sin 2 (Kr B +tf). (12) 

It can be easily proven that the c A (r) and cs(r) fac- 
tors are invariant under sixfold and threefold rotations 
respectively. Let us remember that for intravalley scat- 
tering only, the weights of the density of atomic orbitals 
would be the same on each atoms. 

The isotropic spatial dependence in the LDOS given 
by the Bessel functions in Eqs. (l9llOI) describes long wave- 
length oscillations, since the characteristic wavenumber 
k = e/hvp is much smaller than the cutoff k c . Conse- 
quently, instead of r A and we can use here the con- 
tinous variable r measuring the distance from the impu- 
rity. For distances large enough to satisfy \k\r 3> 1 we 



can use the leading asymptotic expressions Ji(x)Yi(x) — 
cos(2x)/7rx = — Jo(x)Y (x) for x — > oo. Due to this sign 
change on the two sublattices, it is useful to define a com- 
posit short wavelength pattern c(r) = — c A (r) if r is near 
an A site, and c(r) = ce(r) if r is near a B site. Then 
the change of the LDOS due to an impurity at the origin 
can be given by the following compact formula: 



Ap(e,r) = u Q c{r)pl(e) 



cos(2fcr) 
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indicating a long wavelength oscillation with an de- 
cay and a short wavelength spatial pattern given by c(r), 
and shown in Fig. [TJ Although a very rich structure can 
be seen on the plot, the overall threefold rotational sym- 
metry is apparent, and appears to have been observed 
experimentally 15 . 
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FIG. 1: Short wavelength spatial dependence of the LDOS 
due to an impurity at the center. 3D plot of ca(y) + cs(r) = 
|c(r)| is viewed from above, and the color code indicates the 
opposite sign on the two sublattices. 

Due to finite resolution however, the STM tip will mea- 
sure a spatial average of the pattern on Fig. [T] In the 
simpler case of an extended impurity not producing in- 
tervalley scattering, c(r) has atomic densities with equal 
weight on each site of one sublattice and the same weight 
with opposite sign on each site of the other sublattice. 
Clearly, a resolution worse than an elementary cell of 
graphene will lead to cancellation of the leading r _1 de- 
cay in Eq. (fT5|) . leaving us with the next to leading r~ 2 
decay of LDOS. On the other hand for a short range 
impurity potential, internodal scattering contributes as 
well, and yields the short wavelength spatial pattern on 
Fig. [TJ where no cancellation within one unit cell occurs. 
However it is easily shown, that averaging c(r) over three 
neighboring unit cells again leads to cancellation of the 
weights of the atomic densities. Indeed if the weights in 
c(r) are added in the unit cells given by lattice vectors 
R, R + ai and R + ai — &2 (far from the impurity $ 
is the same for all three lattice sites), we obtain zero. 
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Therefore an STM with resolution worse than three ele- 
mentary cells will not be able to measure the leading r 
decay either. It is easily shown that the weights of the 
atomic densities cancel on the six sites of any hexagon of 
the honeycomb lattice. 

The change in the particle density at zero temperature 
due to the impurity is now easily obtained by integrating 
the LDOS up to the Fermi energy £p. Since Ap(e,r) is 
an odd function of energy, An(r) is even in ep, therefore 
we integrate up to — \ep |. As long as kpr ^> 1, where 
kp = \ep\/hvp, we can use the asymptotic forms of the 
Bessel functions to obtain the leading term of the FO as 

. , . UqA c , . , .sm(2kpr) . 
An r = -^e(r) Po ( £F ) / , (14) 
47r r z 

where we disregarded the oscillations due to the cutoff at 
the lower limit of integration as unphysical, since these 
are too fast to be taken seriously in our scheme. The 
FO with wavenumber 2kp decays on both sublattice as 
r~ 2 , but as it was discussed in case of the LDOS, limited 
experimental resolution leads to cancellation of the lead- 
ing power laws, and only the r~ 3 decay will be observed. 
For resolution of about a unit cell this happens for an 
extended impurity not producing internodal scattering, 
but somewhat worse resolution of about three unit cells 
or a hexagon of the honeycomb lattice leads to cancella- 
tion even for a pointlike impurity producing internodal 
scattering as well. Finally it is worth mentioning that 
for half filled graphene, i.e. ep = 0, we obtain r~ 3 decay 
without 2/cf oscillations on both sublattices. 

In conclusion, we have presented a clear and concise 
calculation of the change in the local density of states, 
and the resulting Friedel oscillations in the electron den- 
sity in graphene, due to a single nonmagnetic substitu- 
tional impurity with short range scattering potential. In 
order to achieve atomic scale description we used tight- 
binding wavefunctions with atomic orbitals. As a con- 



sequence of this approach, the short wavelength spatial 
pattern in our results for both LDOS and FO, described 
by c(r) in Eqs. (|13ll4[) , can be considered exact. Using a 
linearized electronic spectrum up to a cutoff affects only 
the rest of the spatial dependence, but for energies close 
to the Dirac point and for distances far from the impu- 
rity, the long wavelength oscillating parts are excellent 
approximations as well. In particular, the LDOS decays 
as r _1 and the FO with wavenumber 2kp as r~ 2 on both 
sublattices. The short wavelength spatial pattern, de- 
picted on Fig. [TJ shows the required rotational symme- 
tries. Since c(r) has alternating signs on neighboring lat- 
tice sites, experimental resolution worse than about three 
unit cells (or a hexagon of the lattice) will lead to cancel- 
lation of the leading power law decays in both quantities 
yielding the next to leading r~ 2 and r~ 3 behavior for 
LDOS and FO respectively. Within the present frame- 
work, the case of an extended defect can also be consid- 
ered by restricting the calculation to intravalley scatter- 
ing by the impurity. This affects only the short wave- 
length spatial pattern c(r), making the weights of the 
atomic densities in it equal (but still of opposite sign on 
the two sublattices) . Consequently, the above mentioned 
cancellation of the leading power law decays occurs al- 
ready at resolution worse than one unit cell. Therefore if 
the experimental resolution falls between one and three 
unit cells, distinction can be made between localized and 
extended defects, since the observed power law decay for 
LDOS for example should follow for the former and 
r~ 2 for the latter. 
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